package fem2.element;

import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.enu.IntegrationType;
import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.math.LinAlg;
import math2.LinAlgUtilities;

/**
 * Linear kinematics structural element using B-bar method
 * 
 * TODO: this element give oscillatory result for strip footing example. TODO
 * need check
 * 
 * @author hbui
 * 
 */
public class BBarStructuralElement extends StructuralElement {

	private IMatrix Bvol = null;

	public BBarStructuralElement(MeshPart mp, MaterialModel mm) {
		super(mp, mm);
	}

	@Override
	public void setIntegrationType(IntegrationType integrationType) {
		if (integrationType != IntegrationType.FULL_INTEGRATION) {
			throw new Error("B-Bar method only support full integration type");
		} else {
			super.setIntegrationType(integrationType);
		}
	}

	/**
	 * Computes the B-Bar operator
	 * 
	 * For 2d:
	 * 
	 * <pre>
	 * TODO: change comment here
	 *     | dN1/dX1  0        dN2/dX1  0       ... |
	 * B = |   0      dN1/dX2    0	    dN2/dX2 ... |
	 *     | dN1/dX2  dN1/dX1  dN2/dX2  dN2/dX1 ... |
	 * </pre>
	 * 
	 * */
	@Override
	public void computeB(MeshPart mp, IMatrix b, double... xi) {
		int n = mp.getBasis().getDimension();
		int dim = mp.getDimension();
		double[][] gGradients = mp.gGradientsAt(xi);
		if (dim == 2) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];

				b.set(0, 2 * i, g1 / 2 + Bvol.get(i, 0));
				b.set(0, 2 * i + 1, -g2 / 2 + Bvol.get(i, 1));
				b.set(1, 2 * i, -g1 / 2 + Bvol.get(i, 0));
				b.set(1, 2 * i + 1, g2 / 2 + Bvol.get(i, 1));
				b.set(2, 2 * i, g2);
				b.set(2, 2 * i + 1, g1);
			}
		} else if (dim == 3) {
			for (int i = 0; i < n; i++) {
				double g1 = gGradients[i][0];
				double g2 = gGradients[i][1];
				double g3 = gGradients[i][2];

				b.set(0, 3 * i, 2 * g1 / 3 + Bvol.get(i, 0));
				b.set(0, 3 * i + 1, -g2 / 3 + Bvol.get(i, 1));
				b.set(0, 3 * i + 2, -g3 / 3 + Bvol.get(i, 2));
				b.set(1, 3 * i, -g1 / 3 + Bvol.get(i, 0));
				b.set(1, 3 * i + 1, 2 * g2 / 3 + Bvol.get(i, 1));
				b.set(1, 3 * i + 2, -g3 / 3 + Bvol.get(i, 2));
				b.set(2, 3 * i, -g1 / 3 + Bvol.get(i, 0));
				b.set(2, 3 * i + 1, -g2 / 3 + Bvol.get(i, 1));
				b.set(2, 3 * i + 2, 2 * g3 / 3 + Bvol.get(i, 2));

				b.set(3, 3 * i, g2);
				b.set(3, 3 * i + 1, g1);
				b.set(4, 3 * i + 1, g3);
				b.set(4, 3 * i + 2, g2);
				b.set(5, 3 * i, g3);
				b.set(5, 3 * i + 2, g1);
			}
		}
	}

	@Override
	public IVector computeStrain(double... xi) {
		IVector e = computeRegularStrain(xi);

		MeshPart mp = getMeshPart();
		int dim = mp.getDimension();
		double ekk = LinAlgUtilities.trace(e);

		double[] uHat = mp.getUHatLocal();
		double omega = LinAlg.dot(uHat.length, uHat, ((Array1DMatrix) Bvol).getValues());

		double tmp = omega - ekk / dim;
		for (int i = 0; i < dim; i++) {
			e.add(i, tmp);
		}

		return e;
	}

	@Override
	public void Initialize() {
		/*
		 * pre-compute Bvol operator
		 */
		MeshPart mp = getMeshPart();
		int dim = mp.getDimension();
		int n = mp.getBasis().getDimension();
		Bvol = Array1DMatrix.createRowMajor(n, dim);
		double t = getMaterialModel().getSectionProperty();
		int ng = getNumIntegrationPoints();
		double Ve = 0.0;

		for (int i = 0; i < ng; i++) {
			double w = getIntegrationWeight(i);
			double[] xi = getIntegrationPoint(i);
			double[][] jac = mp.jacXAt(xi);
			double detJ = LinAlgUtilities.determinantNxN(jac);
			double coeff = w * detJ * t;
			Ve += coeff;

			double[][] gGradients = mp.gGradientsAt(xi);
			for (int j = 0; j < n; j++) {
				for (int k = 0; k < dim; k++) {
					Bvol.add(j, k, coeff * gGradients[j][k]);
				}
			}
		}
		double tmp = 1.0 / (dim * Ve);
		BLAM.scale(tmp, Bvol);
		/*
		 * initialize material and stresses
		 */
		super.Initialize();
	}
}
